#include "IGBT_device_level.h"

void IGBT_device_level::begin()
{
    // IGBT parameters
    VP = 6;
    K = 10;
    KLM = 0;
    A_FET = 0.4;
    M_FET = 1.2;
    N_FET = 1.7;
    BN = 15;
    M_BJT = 2;
    ISAT_BJT = 1E-9;

    // BJT parameters
    k = 1.381E-23;
    q = 1.602E-19;
    T = 300;
    VT = k * T / q;

    // FWD parameters
    M_FWD = 1.21855;
    ISAT_FWD = 1.001e-6;
    RB_FWD = 0.0001;
    INOM = 3000;

    // C_GE
    C_GE = 0.7e-7;
    GCGS = C_GE / h;

    // C_CE
    C0_CE = 2e-8;
    VDIFF_CE = 0.2;
    VSHIFT_CE = 0;
    ALPHA_CE = 0.5;
    BETA_CE = 2;
    DELTA_CE = 0.0001;

    // C_CG
    C0_CG = 5e-8;
    VDIFF_CG_OFF = 0.1;
    VDIFF_CG_ON = 0.1;
    VSHIFT_CG_OFF = 0;
    VSHIFT_CG_ON = 0;
    ALPHA_CG_OFF = 0.5;
    ALPHA_CG_ON = 0.5;
    BETA_CG_OFF = 2;
    BETA_CG_ON = 2;
    DELTA_CG_OFF = 0.0001;
    DELTA_CG_ON = 0.0001;

    // Ceb(diff)
    TAU_BE = 5e-8;

    // connector
    Rg = 10.2;
    LG = 0;
    LC = 5e-9;
    RC = 1.15625e-5;
    LE = 5e-9;
    RE = 1.15625e-5;
    RAUX = 1.15625e-5;
    LAUX = 1.75e-10;
}

void IGBT_device_level::dsdyn_in_iteration(double delt, double simulation_time)
{
    this->delt = delt;
    this->h = this->delt;
    this->simulation_time = simulation_time;

    if (simple)
    {
        if(Vs > 0)
        {
            Gout = 1.0e3;
            Jout = 0;
        }
        else
        {
            Gout = 1.0e-6;
            Jout = 0;
        }
        Itotal = -2333;
        return;
    }

    // Do not delete the comments in this function！

    // double VP, K, KLM, A_FET, M_FET, N_FET, BN, M_BJT, ISAT_BJT, Vsat, Isat, VDS, ID;
    // double GmosVGS, GmosVDS, Imoseq, ID_Factor1, ID_Factor2, ID_Factor3, d1VDS, d2VDS, d3VDS, d2VGS, d3VGS;
    // double k, q, T, VT, VEB, IB, Gbjt, IBeq, RP_BJT, M_FWD, ISAT_FWD, IFWD, GFWD, IFWDeq, VFVD;
    // double C_GE, GCGS, ICGS_eq, Rg;
    // VP = 6;
    // K = 10;
    // KLM = 0;
    // A_FET = 0.4;
    // M_FET = 1.2;
    // N_FET = 1.7;
    // BN = 15;
    // M_BJT = 2;
    // ISAT_BJT = 1E-9;

    VGS = V_nodeG - V_nodeS;

    // 2.2 MOSFET
    Vsat = A_FET * pow((VGS - VP), M_FET);
    Isat = K / 2 * pow((VGS - VP), N_FET);
    VDS = V_nodeD - V_nodeS;

    if (VDS < 0)
    {
        VDS = 0;
    }

    if (VGS <= VP)
    {
        ID = 0;

        GmosVGS = 0;
        GmosVDS = 0;
        Imoseq = 0;
    }
    else if (VDS < Vsat)
    {
        ID_Factor1 = 1 + KLM * VDS;
        ID_Factor2 = 2 - VDS / Vsat;
        ID_Factor3 = VDS / Vsat;
        ID = Isat * ID_Factor1 * ID_Factor2 * ID_Factor3;

        d1VDS = KLM;
        d2VDS = -1 / Vsat;
        d3VDS = 1 / Vsat;
        GmosVDS = Isat * (d1VDS * ID_Factor2 * ID_Factor3 + ID_Factor1 * d2VDS * ID_Factor3 + ID_Factor1 * ID_Factor2 * d3VDS);
        d2VGS = VDS / pow(Vsat, 2) * A_FET * M_FET * pow((VGS - VP), (M_FET - 1));
        d3VGS = -VDS / pow(Vsat, 2) * A_FET * M_FET * pow((VGS - VP), (M_FET - 1));
        GmosVGS = K / 2 * N_FET * pow((VGS - VP), (N_FET - 1)) * ID_Factor1 * ID_Factor2 * ID_Factor3 + Isat * ID_Factor1 * d2VGS * ID_Factor3 + Isat * ID_Factor1 * ID_Factor2 * d3VGS;
        Imoseq = ID - GmosVDS * VDS - GmosVGS * (VGS);
    }
    else
    {
        ID = Isat * (1 + KLM * VDS);

        GmosVDS = Isat * KLM;
        GmosVGS = K / 2 * N_FET * pow((VGS - VP), (N_FET - 1)) * (1 + KLM * VDS);
        Imoseq = ID - GmosVDS * VDS - GmosVGS * (VGS);
    }

    // 2.3 BJT
    // k = 1.381E-23;
    // q = 1.602E-19;
    // T = 300;
    // VT = k * T / q;
    VEB = V_node1 - V_nodeD;

    if (VEB > 2)
    {
        VEB = 2;
    }

    IB = ISAT_BJT * (exp(VEB / M_BJT / VT) - 1);
    Gbjt = ISAT_BJT * exp(VEB / M_BJT / VT) / M_BJT / VT;
    IBeq = IB - Gbjt * VEB;

    RP_BJT = 125;

    ////////////2.4 FWD

    // M_FWD = 1.21855;
    // ISAT_FWD = 1.001e-6;
    VFVD = V_nodeAUX - V_nodemidFWD;
    ///////////////////////////////////important////////////////////////////////////
    if (VFVD > 1)
    {
        VFVD = 1;
    }
    ////////////////////////////////////////////////////////////////////////////////
    IFWD = ISAT_FWD * (exp(VFVD / M_FWD / VT) - 1);
    GFWD = ISAT_FWD * exp(VFVD / M_FWD / VT) / M_FWD / VT;
    IFWDeq = IFWD - GFWD * VFVD;

    // double RB_FWD = 0.0001;
    // double INOM = 3000;
    double R_FWD = RB_FWD / sqrt(1 + IFWD / INOM);
    double G_FWD = 1 / R_FWD;

    /////////////2.5 C_GE
    // C_GE = 0.7e-7;
    // GCGS = C_GE / h;
    ICGS_eq = -GCGS * VGS_old;

    /////////////2.6 C_CE
    // double C0_CE = 2e-8;
    // double VDIFF_CE = 0.2;
    // double VSHIFT_CE = 0;
    // double ALPHA_CE = 0.5;
    // double BETA_CE = 2;
    // double DELTA_CE = 0.0001;

    // double VCE, B, VVJNCT_CE, CEa, CEb, AA_CE, BB_CE;

    VCE = V_node1 - V_nodeS; // QCE_old
    B = VSHIFT_CE * VDIFF_CE;
    VVJNCT_CE = B - VCE;
    CEa = BETA_CE - 1;
    CEb = 1 - DELTA_CE;

    AA_CE = (ALPHA_CE * CEb / CEa) / VDIFF_CE;
    BB_CE = (1 - DELTA_CE) * VDIFF_CE / (1 - ALPHA_CE);

    // double CCE, ICE, GCE, ICEeq;
    if (VVJNCT_CE >= 0)
    {
        CCE = C0_CE * (1 + CEa * (1 - exp(-VVJNCT_CE * AA_CE)));

        // Charge Calculation by Integration of CCE(from 0 to VCE)
        QCE = C0_CE * (VCE + CEa * (VCE - exp(VCE * AA_CE) / AA_CE)) + C0_CE * CEa / AA_CE;
        ICE = (QCE - QCE_old) / h;
        GCE = CCE / h;
        ICEeq = ICE - GCE * VCE;
    }

    else
    {
        CCE = C0_CE * (DELTA_CE + CEb / pow((1 - VVJNCT_CE / VDIFF_CE), ALPHA_CE));
        // Charge Calculation by Integration of CCE(from 0 to VCE)
        QCE = C0_CE * (DELTA_CE * VCE + BB_CE * pow((1 + VCE / VDIFF_CE), (1 - ALPHA_CE))) - C0_CE * BB_CE;
        ICE = (QCE - QCE_old) / h;
        GCE = CCE / h;
        ICEeq = ICE - GCE * VCE;
    }
    ///////////////////2.7 C_CG
    // double C0_CG = 5e-8;
    // double VDIFF_CG_OFF = 0.1;
    // double VDIFF_CG_ON = 0.1;
    // double VSHIFT_CG_OFF = 0;
    // double VSHIFT_CG_ON = 0;
    // double ALPHA_CG_OFF = 0.5;
    // double ALPHA_CG_ON = 0.5;
    // double BETA_CG_OFF = 2;
    // double BETA_CG_ON = 2;
    // double DELTA_CG_OFF = 0.0001;
    // double DELTA_CG_ON = 0.0001;

    // double VCG, VVJNCT_CG, CGa, CGb, AA_CG, BB_CG, CCG, ICG, GCG, ICGeq;

    VCG = V_node1 - V_nodeG; // QCG_old
    if (VGS <= VP)
    {
        B = VSHIFT_CG_OFF * VDIFF_CG_OFF;
        VVJNCT_CG = B - VCG;
        CGa = BETA_CG_OFF - 1;
        CGb = 1 - DELTA_CG_OFF;

        AA_CG = ALPHA_CG_OFF * CGb / CGa / VDIFF_CG_OFF;
        BB_CG = (1 - DELTA_CG_OFF) * VDIFF_CG_OFF / (1 - ALPHA_CG_OFF);

        if (VVJNCT_CG >= 0)
        {
            CCG = C0_CG * (1 + CGa * (1 - exp(-VVJNCT_CG * AA_CG)));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (VCG + CGa * (VCG - exp(VCG * AA_CG) / AA_CG)) + C0_CG * CGa / AA_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
        else
        {
            CCG = C0_CG * (DELTA_CG_OFF + CGb / pow((1 - VVJNCT_CG / VDIFF_CG_OFF), ALPHA_CG_OFF));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (DELTA_CG_OFF * VCG + BB_CG * pow((1 + VCG / VDIFF_CG_OFF), (1 - ALPHA_CG_OFF))) - C0_CG * BB_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
    }
    else
    {
        B = VSHIFT_CG_ON * VDIFF_CG_ON;
        VVJNCT_CG = B - VCG;
        CGa = BETA_CG_ON - 1;
        CGb = 1 - DELTA_CG_ON;

        AA_CG = ALPHA_CG_ON * CGb / CGa / VDIFF_CG_ON;
        BB_CG = (1 - DELTA_CG_ON) * VDIFF_CG_ON / (1 - ALPHA_CG_ON);

        if (VVJNCT_CG >= 0)
        {
            CCG = C0_CG * (1 + CGa * (1 - exp(-VVJNCT_CG * AA_CG)));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (VCG + CGa * (VCG - exp(VCG * AA_CG) / AA_CG)) + C0_CG * CGa / AA_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
        else
        {
            CCG = C0_CG * (DELTA_CG_ON + CGb / pow((1 - VVJNCT_CG / VDIFF_CG_ON), ALPHA_CG_ON));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (DELTA_CG_ON * VCG + BB_CG * pow((1 + VCG / VDIFF_CG_ON), (1 - ALPHA_CG_ON))) - C0_CG * BB_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
    }

    //////////////2.8 Ceb(diff)
    // double TAU_BE = 5e-8;
    // double Ceb, Ieb, GCeb, ICebeq;

    Ceb = TAU_BE * (IB + ISAT_BJT) / M_BJT / VT;
    Qeb = TAU_BE * ISAT_BJT * (exp(VEB / M_BJT / VT) - 1);
    Ieb = (Qeb - Qeb_old) / h;
    GCeb = Ceb / h;
    ICebeq = Ieb - GCeb * VEB;

    /////////////2.6 Connector

    // Rg = 10.2;
    // double LG = 0;
    // double LC = 5e-9;
    // double RC = 1.15625e-5;
    // double LE = 5e-9;
    // double RE = 1.15625e-5;
    // double RAUX = 1.15625e-5;
    // double LAUX = 1.75e-10;

    f[0] = 0;                                         // node C
    f[1] = -(BN + 1) * IBeq - ICEeq - ICGeq - ICebeq; // node 1
    f[2] = -Imoseq + IBeq + ICebeq;                   // node D
    f[3] = Vs / Rg - ICGS_eq + ICGeq;                 // node G
    f[4] = Imoseq + BN * IBeq + ICGS_eq + ICEeq;      // node S
    f[5] = -IFWDeq - Vs / Rg;                         // node AUX
    f[6] = IFWDeq;                                    // node mid FWD
    f[7] = 0;                                         // node E

    G[0][0] = 1 / RC;
    G[0][1] = -1 / RC;
    G[0][2] = 0;
    G[0][3] = 0;
    G[0][4] = 0;
    G[0][5] = 0;
    G[0][6] = 0;
    G[0][7] = 0;

    G[1][0] = -1 / RC;
    G[1][1] = (BN + 1) * Gbjt + 1 / RP_BJT + G_FWD + 1 / RC + GCE + GCG + GCeb;
    G[1][2] = -(BN + 1) * Gbjt - 1 / RP_BJT - GCeb;
    G[1][3] = -GCG;
    G[1][4] = -GCE;
    G[1][5] = 0;
    G[1][6] = -G_FWD;
    G[1][7] = 0;

    G[2][0] = 0;
    G[2][1] = -Gbjt - 1 / RP_BJT - GCeb;
    G[2][2] = GmosVDS + Gbjt + 1 / RP_BJT + GCeb;
    G[2][3] = GmosVGS;
    G[2][4] = -GmosVDS - GmosVGS;
    G[2][5] = 0;
    G[2][6] = 0;
    G[2][7] = 0;

    G[3][0] = 0;
    G[3][1] = -GCG;
    G[3][2] = 0;
    G[3][3] = 1 / Rg + GCGS + GCG;
    G[3][4] = -GCGS;
    G[3][5] = -1 / Rg;
    G[3][6] = 0;
    G[3][7] = 0;

    G[4][0] = 0;
    G[4][1] = -BN * Gbjt - GCE;
    G[4][2] = -GmosVDS + BN * Gbjt;
    G[4][3] = -GmosVGS - GCGS;
    G[4][4] = GmosVDS + GmosVGS + GCGS + 1 / RAUX + GCE;
    G[4][5] = -1 / RAUX;
    G[4][6] = 0;
    G[4][7] = 0;

    G[5][0] = 0;
    G[5][1] = 0;
    G[5][2] = 0;
    G[5][3] = -1 / Rg;
    G[5][4] = -1 / RAUX;
    G[5][5] = GFWD + 1 / Rg + 1 / RAUX + 1 / RE;
    G[5][6] = -GFWD;
    G[5][7] = -1 / RE;

    G[6][0] = 0;
    G[6][1] = -G_FWD;
    G[6][2] = 0;
    G[6][3] = 0;
    G[6][4] = 0;
    G[6][5] = -GFWD;
    G[6][6] = G_FWD + GFWD;
    G[6][7] = 0;

    G[7][0] = 0;
    G[7][1] = 0;
    G[7][2] = 0;
    G[7][3] = 0;
    G[7][4] = 0;
    G[7][5] = -1 / RE;
    G[7][6] = 0;
    G[7][7] = 1 / RE;

    for (int i = 0; i < 8; i++)
    {
        J8(i) = f[i];
        for (int j = 0; j < 8; j++)
        {
            G8(i, j) = G[i][j];
        }
    }

    Y11 = G8.topLeftCorner<1, 1>();
    Y12 = G8.block<1, 6>(0, 1);
    Y21 = Y12.transpose();
    Y22 = G8.block<6, 6>(1, 1);
    Y22_inv = Y22.inverse();

    Jex(0) = J8(0);
    Jin = J8.block<6, 1>(1, 0);

    Geq = Y11 - Y12 * Y22_inv * Y21;
    Jeq = Jex - Y12 * Y22_inv * Jin;

    Gout = Geq(0);
    Jout = Jeq(0) * 1.0e-3; // 外部电流单位为kA，内部电流单位为A

    Itotal = (V_nodeC - V_node1) / RC;
    return;
}

void IGBT_device_level::dsout_in_iteration()
{
    Vin = Y22_inv * (Jin - Y21 * Vex);
    V_IGBT[0] = Vex(0);
    for (int i = 1; i < 7; i++)
    {
        V_IGBT[i] = Vin(i - 1);
    }
    V_IGBT[7] = 0;
}

void IGBT_device_level::dsout_out_iteration()
{
    VGS_old = VGS;
    QCE_old = QCE;
    QCG_old = QCG;
    Qeb_old = Qeb;
}